Wykorzystane biblioteki:

library(ggplot2) #Do tworzenia regularnych wykresów
library(plotly) #Do tworzenia interaktywnych wykresów
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(dplyr) #Do manipulacji na danych
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(magrittr) #Dla operatora %>%
library(caTools)
library(reshape2) #Dla funkcji 'melt' upraszczającej stworzenie heatmapy korelacji
library(knitr) #Dla prezentacji dataframe'a w sensowny sposób w RMarkdown

Wymuszanie identycznych rezultatów dla kolejnych wywołań i wczytanie danych:

set.seed(2112) #Zadbanie o taki sam rezultat dla kolejnych wywołań
data <- read.csv('Life_Expectancy_Data.csv')

Funkcje pomocniczne do pokazania podsumowań:

sum_nas <- function(x, na.rm){
    sum(is.na(x))
}
statistics <- list(list('Nothing_to_see_here', "Column name"), 
                   list(mean, "Mean"), list(sd, "Standard Deviation"),
                   list(min, "Minimum"), list(max, "Maximum"),
                  list(median, "Median"), list(sum_nas, "Amount of NaN"))

all_cols = list()
for (x in statistics){
    all_cols[[length(all_cols)+1]] <- x[[2]]
}

df_summary <- data.frame(matrix(ncol = length(all_cols), nrow = 0))
colnames(df_summary) <- all_cols

Przedstawienie danych ogólnych o datasecie i zmiennych kategorycznych:

cols <- ncol(data)
cat("Number of rows: ", nrow(data), "\nNumber of columns: ", cols, "\n")
## Number of rows:  2938 
## Number of columns:  22
ite <- 1

for (x in data){
    
    dist = "   "
    if (class(x) != 'factor'){
        new_row <- list()
        for (y in statistics){
            if (y[[2]] == 'Column name')
                new_row[[length(new_row)+1]] <- colnames(data)[ite]
            else
                new_row[[length(new_row)+1]] <- y[[1]](x, na.rm=TRUE)
        }
        df_summary[nrow(df_summary) + 1, ] <- new_row
    }
    else{
        names <- unique(x)
        cat('\nFactor column: ', colnames(data)[ite], ', Unique values: ', paste(shQuote(names, type="cmd"), collapse=", "), '\n')
    }
    ite <- ite+1
}
## 
## Factor column:  Country , Unique values:  "Afghanistan", "Albania", "Algeria", "Angola", "Antigua and Barbuda", "Argentina", "Armenia", "Australia", "Austria", "Azerbaijan", "Bahamas", "Bahrain", "Bangladesh", "Barbados", "Belarus", "Belgium", "Belize", "Benin", "Bhutan", "Bolivia (Plurinational State of)", "Bosnia and Herzegovina", "Botswana", "Brazil", "Brunei Darussalam", "Bulgaria", "Burkina Faso", "Burundi", "Cote d'Ivoire", "Cabo Verde", "Cambodia", "Cameroon", "Canada", "Central African Republic", "Chad", "Chile", "China", "Colombia", "Comoros", "Congo", "Cook Islands", "Costa Rica", "Croatia", "Cuba", "Cyprus", "Czechia", "Democratic People's Republic of Korea", "Democratic Republic of the Congo", "Denmark", "Djibouti", "Dominica", "Dominican Republic", "Ecuador", "Egypt", "El Salvador", "Equatorial Guinea", "Eritrea", "Estonia", "Ethiopia", "Fiji", "Finland", "France", "Gabon", "Gambia", "Georgia", "Germany", "Ghana", "Greece", "Grenada", "Guatemala", "Guinea", "Guinea-Bissau", "Guyana", "Haiti", "Honduras", "Hungary", "Iceland", "India", "Indonesia", "Iran (Islamic Republic of)", "Iraq", "Ireland", "Israel", "Italy", "Jamaica", "Japan", "Jordan", "Kazakhstan", "Kenya", "Kiribati", "Kuwait", "Kyrgyzstan", "Lao People's Democratic Republic", "Latvia", "Lebanon", "Lesotho", "Liberia", "Libya", "Lithuania", "Luxembourg", "Madagascar", "Malawi", "Malaysia", "Maldives", "Mali", "Malta", "Marshall Islands", "Mauritania", "Mauritius", "Mexico", "Micronesia (Federated States of)", "Monaco", "Mongolia", "Montenegro", "Morocco", "Mozambique", "Myanmar", "Namibia", "Nauru", "Nepal", "Netherlands", "New Zealand", "Nicaragua", "Niger", "Nigeria", "Niue", "Norway", "Oman", "Pakistan", "Palau", "Panama", "Papua New Guinea", "Paraguay", "Peru", "Philippines", "Poland", "Portugal", "Qatar", "Republic of Korea", "Republic of Moldova", "Romania", "Russian Federation", "Rwanda", "Saint Kitts and Nevis", "Saint Lucia", "Saint Vincent and the Grenadines", "Samoa", "San Marino", "Sao Tome and Principe", "Saudi Arabia", "Senegal", "Serbia", "Seychelles", "Sierra Leone", "Singapore", "Slovakia", "Slovenia", "Solomon Islands", "Somalia", "South Africa", "South Sudan", "Spain", "Sri Lanka", "Sudan", "Suriname", "Swaziland", "Sweden", "Switzerland", "Syrian Arab Republic", "Tajikistan", "Thailand", "The former Yugoslav republic of Macedonia", "Timor-Leste", "Togo", "Tonga", "Trinidad and Tobago", "Tunisia", "Turkey", "Turkmenistan", "Tuvalu", "Uganda", "Ukraine", "United Arab Emirates", "United Kingdom of Great Britain and Northern Ireland", "United Republic of Tanzania", "United States of America", "Uruguay", "Uzbekistan", "Vanuatu", "Venezuela (Bolivarian Republic of)", "Viet Nam", "Yemen", "Zambia", "Zimbabwe" 
## 
## Factor column:  Status , Unique values:  "Developing", "Developed"

Przedstawienie informacji o zmiennych numerycznych dataseta:

#Wzorzec: https://www.rdocumentation.org/packages/knitr/versions/1.36/topics/knit_print
knit_print.data.frame = function(x, ...) {
    res = paste(c("", "", kable(x, output = FALSE)), collapse = "\n")
    asis_output(res)
}
knit_print(df_summary)
Column name Mean Standard Deviation Minimum Maximum Median Amount of NaN
Year 2.007519e+03 4.613841e+00 2000.00000 2.015000e+03 2.008000e+03 0
Life.expectancy 6.922493e+01 9.523867e+00 36.30000 8.900000e+01 7.210000e+01 10
Adult.Mortality 1.647964e+02 1.242921e+02 1.00000 7.230000e+02 1.440000e+02 10
infant.deaths 3.030395e+01 1.179265e+02 0.00000 1.800000e+03 3.000000e+00 0
Alcohol 4.602861e+00 4.052413e+00 0.01000 1.787000e+01 3.755000e+00 194
percentage.expenditure 7.382513e+02 1.987915e+03 0.00000 1.947991e+04 6.491291e+01 0
Hepatitis.B 8.094046e+01 2.507002e+01 1.00000 9.900000e+01 9.200000e+01 553
Measles 2.419592e+03 1.146727e+04 0.00000 2.121830e+05 1.700000e+01 0
BMI 3.832125e+01 2.004403e+01 1.00000 8.730000e+01 4.350000e+01 34
under.five.deaths 4.203574e+01 1.604455e+02 0.00000 2.500000e+03 4.000000e+00 0
Polio 8.255019e+01 2.342805e+01 3.00000 9.900000e+01 9.300000e+01 19
Total.expenditure 5.938190e+00 2.498320e+00 0.37000 1.760000e+01 5.755000e+00 226
Diphtheria 8.232408e+01 2.371691e+01 2.00000 9.900000e+01 9.300000e+01 19
HIV.AIDS 1.742104e+00 5.077784e+00 0.10000 5.060000e+01 1.000000e-01 0
GDP 7.483158e+03 1.427017e+04 1.68135 1.191727e+05 1.766948e+03 448
Population 1.275338e+07 6.101210e+07 34.00000 1.293859e+09 1.386542e+06 652
thinness..1.19.years 4.839704e+00 4.420195e+00 0.10000 2.770000e+01 3.300000e+00 34
thinness.5.9.years 4.870317e+00 4.508882e+00 0.10000 2.860000e+01 3.300000e+00 34
Income.composition.of.resources 6.275511e-01 2.109036e-01 0.00000 9.480000e-01 6.770000e-01 167
Schooling 1.199279e+01 3.358920e+00 0.00000 2.070000e+01 1.230000e+01 163

Zadbanie o to, aby wykresy były niesamowicie eleganckie:

options(repr.plot.width=24, repr.plot.height=16)

Tworzenie eleganckich wykresów rozkładu zmiennych: warto zwrócić uwagę na wykorzystanie density-plotu dla danych ciągłych i barplota dla danych dyskretnych. Uwzględniono - dla kompletności - rozkład krajów i lat.

for (x in colnames(data)){
    if (is.numeric(data[[x]]) && x!='Year')
        print(data %>%
              filter(!is.na(data[[x]])) %>%
              ggplot(aes_string(x=x)) + geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8))
    else
        print(data %>%
              filter(!is.na(data[[x]])) %>%
              ggplot(aes_string(x=x)) + geom_bar() + 
                  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 6, hjust = 1)))
}

Dodatkowe wykresy prezentujący korelacje pomiędzy Life.expectancy a pozostałymi zmiennymi: pozwalają one zauważyć, że pewne wartości mogły nie zostać poprawnie udokumentowane.

for (x in colnames(data)){
    if (x == 'Life.expectancy')
        next
    plot_data <- data %>% filter(!is.na(data[[x]]) & !is.na(data[['Life.expectancy']]))
    plot_itself <- plot_data %>% ggplot(aes_string(x=x, y='Life.expectancy')) + geom_point()
    if (x == 'Country')
        plot_itself <- plot_itself + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 6, hjust = 1))
    print(plot_itself)
}


Z wykresów tych można wywnioskować, że część obsserwacji używała innej skali dla pewnych zmiennych niż wymagana: np. BMI rzędu poniżej 7 (które występowało w krajach o relatywnie wysokiej długości życia) oznacza wagę rzędu ~10 kg dla człowieka o wzroście 1.70 m - można powątpiewać w zasadność tych liczb. Podobnie ma się rzecz z wyszczepialnością; W kolumnie Diphtheria jest bardzo widoczna linia na wysokości ok. x=10; istnieją co najmniej 2 kraje, dla których wartość ‘Diphtheria’ “skakała” pomiędzy wartościami z przedziału 8-9 i 80-90 - Norwegia i Mołdawia.

knit_print(data %>% filter(Country == 'Republic of Moldova' | Country == 'Norway') %>% select(c(Year, Country, Diphtheria, BMI)) %>% arrange(Country, Year))
Year Country Diphtheria BMI
2000 Norway 9 53.3
2001 Norway 91 54.0
2002 Norway 93 54.6
2003 Norway 92 55.3
2004 Norway 92 55.9
2005 Norway 91 56.5
2006 Norway 94 57.0
2007 Norway 93 57.5
2008 Norway 94 58.0
2009 Norway 94 58.5
2010 Norway 93 58.9
2011 Norway 94 59.4
2012 Norway 95 59.8
2013 Norway 94 6.3
2014 Norway 93 6.8
2015 Norway 95 61.2
2000 Republic of Moldova 95 46.5
2001 Republic of Moldova 97 46.8
2002 Republic of Moldova 97 47.2
2003 Republic of Moldova 98 47.6
2004 Republic of Moldova 98 47.9
2005 Republic of Moldova 98 48.3
2006 Republic of Moldova 97 48.7
2007 Republic of Moldova 92 49.1
2008 Republic of Moldova 9 49.5
2009 Republic of Moldova 85 49.9
2010 Republic of Moldova 9 5.4
2011 Republic of Moldova 93 5.9
2012 Republic of Moldova 92 51.5
2013 Republic of Moldova 9 52.1
2014 Republic of Moldova 9 52.7
2015 Republic of Moldova 87 53.4

Można podejrzewać, że część wyników podano w błędnych skalach - wiedza ta zostanie wykorzystana przy tworzeniu regresora. W przypadku tych dwóch zmiennych można relatywnie łatwo wyodrębnić dane w innej skali; z kolei dla np. Adult.mortality nie dokonano zmian, gdyż linią podziału nie jest stała.

Heatmapa przedstawiająca korelacje między zmiennymi:

cor_matrix <- data %>% select(-c(Country, Status)) %>% cor(use = "complete.obs") %>% round(2)
melt_matrix <- melt(cor_matrix, na.rm = TRUE)

#Wzorzec I: https://r-charts.com/correlation/heat-map-ggplot2/
#Wzorzec II: http://sthda.com/english/wiki/ggplot2-quick-correlation-matrix-heatmap-r-software-and-data-visualization

ggplot(data = melt_matrix, aes(Var2, Var1, fill = value))+
    geom_tile()+
    scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
        midpoint = 0, limit = c(-1,1), space = "Lab", 
        name="Pearson\nCorrelation") +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, 
        size = 5, hjust = 1))+
    geom_text(aes(label = value), color = 'black', size = 1.5)+
    coord_fixed()

Preparacja danych do stworzenia interaktywnego wykresu długości życia w funkcji kraju i roku:

unique_countries <- unique(data$Country)
buttons_countries <- list()

for (i in seq_along(unique_countries)){
    buttons_countries[[length(buttons_countries)+1]] <- list(method = "restyle",
               args = list("transforms[0].value", unique_countries[i]),
               label = unique_countries[i])
}

Wykres długości życia w funkcji kraju i roku:

#Wzorzec: https://stackoverflow.com/questions/63906441/plotly-r-using-colour-and-transformation-with-a-line-plot
p <- data %>% filter(!is.na(Life.expectancy)) %>%
  plot_ly(
    type = 'scatter', 
    x = ~Year, 
    y = ~Life.expectancy,
    mode = 'markers',
    transforms = list(
      list(
        type = 'filter',
        target = ~Country,
        operation = '=',
        value = unique_countries[1]
      )
  )) %>% layout(
    updatemenus = list(
      list(
        type = 'dropdown',
        active = 0,
        buttons = buttons_countries
      )
    )
  )

p

Preparacja danych dla regresora - w szczególności wykorzystano informację o tym, czy w danej kolumnie zaszedł NaN - tworzono wtedy nową koumnę ‘col_name_is_NaN’ sprawdzającą, czy dana obserwacja ma NaN-a w danej kolumnie, a następnie przypisując medianę istniejących wartości z danej kolumny do NaN-a. Co może się wydawać zaskakujące, to osiągało lepszą skuteczność w regresorze niż wstawienie 0 zamiast NaN-a - może to wynikać z tego, że medianę liczono przed podziałem na zbiór treningowy i testowy. Ponadto wykorzystano obserwację z analizy korelacji między Life.expectancy a pewnymi problematycznie przeskalowanymi zmiennymi.

change <- function(data, name, borderline, multiplier) {
    data[[name]][data[[name]] < borderline & !is.na(data[[name]])] <- data[[name]][data[[name]] < borderline & !is.na(data[[name]])]*multiplier
    data
}

data <- data %>% filter(!is.na(Life.expectancy))

data <- change(data, 'BMI', 8.5, 10)
data <- change(data, 'Diphtheria', 10, 10)
data <- change(data, 'Polio', 10, 10)

for (x in colnames(data)){
    if (data %>% filter(is.na(data[[x]])) %>% nrow == 0)
        next
    data[[paste(x, '_na', sep='')]] <- rep(0, nrow(data))
    data[[paste(x, '_na', sep='')]][is.na(data[[x]])]  <- 1
    mid = data[[x]]
    data[[x]][is.na(data[[x]])]  <- median(data[[x]], na.rm=TRUE)
}

Usunięcie powtarzających się kolumn po dodaniu kolumn z NaN-em:

marked = 'Q'
while (marked != '-'){
    marked='-'
    for (x in colnames(data)){
        for (y in colnames(data)){
            if (!is.numeric(data[[x]]) || !is.numeric(data[[y]]) || x==y)
                next
            summa = ifelse (data[[x]] == data[[y]], 1, 0)
            if (sum(summa) == nrow(data)){
                marked = y
                break
            }
        }
        if (marked!='-') break
    }
    if (marked!='-')
        data <- data %>% select(-c(y))
}
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(y)` instead of `y` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

Wykonanie regresora (warto zauważyć, że regresor nie wykorzystuje kolumny “Country”):

tmp_data= sample.split(data, SplitRatio = 0.3)
train = subset(data,tmp_data==TRUE)
test = subset(data,tmp_data==FALSE)

lm_death <- train %>% select(-c(Country)) %>% lm(formula=Life.expectancy ~ .)
summary(lm_death)
## 
## Call:
## lm(formula = Life.expectancy ~ ., data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.8979  -2.1109   0.0173   2.2059  11.8731 
## 
## Coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         7.970e+01  6.681e+01   1.193  0.23324    
## Year                               -1.487e-02  3.333e-02  -0.446  0.65548    
## StatusDeveloping                   -2.127e+00  4.586e-01  -4.638 4.07e-06 ***
## Adult.Mortality                    -1.780e-02  1.417e-03 -12.564  < 2e-16 ***
## infant.deaths                       7.543e-02  1.568e-02   4.811 1.78e-06 ***
## Alcohol                             5.049e-02  4.438e-02   1.138  0.25564    
## percentage.expenditure              2.256e-04  1.457e-04   1.548  0.12195    
## Hepatitis.B                        -2.800e-02  6.520e-03  -4.294 1.95e-05 ***
## Measles                            -5.050e-07  1.520e-05  -0.033  0.97351    
## BMI                                 9.851e-02  1.330e-02   7.405 3.16e-13 ***
## under.five.deaths                  -5.596e-02  1.139e-02  -4.913 1.08e-06 ***
## Polio                               1.026e-01  3.251e-02   3.156  0.00166 ** 
## Total.expenditure                  -9.412e-03  5.698e-02  -0.165  0.86884    
## Diphtheria                          4.490e-02  3.236e-02   1.387  0.16568    
## HIV.AIDS                           -4.335e-01  3.166e-02 -13.692  < 2e-16 ***
## GDP                                 5.362e-06  2.260e-05   0.237  0.81251    
## Population                         -2.231e-09  2.523e-09  -0.884  0.37698    
## thinness..1.19.years                2.399e-03  7.893e-02   0.030  0.97576    
## thinness.5.9.years                  7.808e-02  7.720e-02   1.011  0.31214    
## Income.composition.of.resources     5.386e+00  1.021e+00   5.273 1.70e-07 ***
## Schooling                           4.919e-01  8.192e-02   6.005 2.83e-09 ***
## Alcohol_na                          2.760e+00  1.420e+00   1.944  0.05225 .  
## Hepatitis.B_na                      7.275e-01  4.204e-01   1.730  0.08391 .  
## BMI_na                             -3.485e+00  1.307e+00  -2.666  0.00783 ** 
## Polio_na                            2.835e+00  2.085e+00   1.360  0.17423    
## Total.expenditure_na               -1.277e+00  1.295e+00  -0.986  0.32453    
## GDP_na                             -9.976e-01  6.503e-01  -1.534  0.12537    
## Population_na                       1.261e+00  5.424e-01   2.325  0.02030 *  
## Income.composition.of.resources_na -1.116e+00  7.520e-01  -1.483  0.13836    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.661 on 848 degrees of freedom
## Multiple R-squared:  0.8592, Adjusted R-squared:  0.8546 
## F-statistic: 184.8 on 28 and 848 DF,  p-value: < 2.2e-16
pred <- predict(lm_death, test)

Można zauważyć, że - z pominięciem kraju - najlepszymi predyktorami zmiennej decyzyjnej były zmienne: Schooling, Hepatitis.B, Status, Adult.Mortality, Diphtheria, HIV.AIDS i Income.composition.of.resources - są to też zmienne, które zostały wskazane jako dobrze skorelowane ze zmienną decyzyjną w pokazanej heatmapie. Zapewne nie można argumentować o wpływie tych zmiennych na oczekiwaną długość życia, ponieważ korelacja może informować jedynie o związku, ale nie o przyczynowości.

Wartość \(rmse\):

sqrt(mean((pred - test$Life.expectancy)^2))
## [1] 3.846439

Po przeskalowaniu problematycznych wartości BMI i części szczepień wartości \(rmse\) spadły z 4.11 do 3.8. Przy wykorzystaniu kraju jako predyktora spadłyby jeszcze bardziej (poniżej 2), natomiast jego wykorzystanie można potraktować jako problematyczne; w praktyce można traktować kraj jako pewne ID i podzielić dane kraj po kraju do danych treningowych albo testowych - wtedy predyktor bardzo silnie zależny od kraju kompletnie by sobie nie radził.